1. Field of the Invention
The present invention relates to a process simulation method using a computer system applicable to semiconductor device fabrication and more particularly, to a computer simulation method of oxidation of silicon (Si), in which the diffusion equation of oxidant is solved to find the surface concentration of the oxidant at the interface between a Si region and a silicon dioxide (SiO.sub.2) region.
2. Description of the Prior Art
A process simulator is a computer system to simulate various processes in the semiconductor device fabrication, such as oxidation, diffusion, etching, and ion-implantation, thereby predicting the details of the resulting device structure, such as profiles of doped impurities and topography of conductive or dielectric materials. If the device structure of a Large-Scale Integrated circuit (LSI) is optimized by the use of the process simulator in such a way that the LSI exhibits the desired electrical characteristics, the developmental cost and period for the LSI can be drastically reduced compared with the case where the LSI is actually fabricated for the purpose of experiments.
Conventionally, the process simulator designed for semiconductor device fabrication is equipped with built-in simulation models applicable to the individual fabrication processes. For example, a simulation model of the time-dependent thickness of an oxide region is disclosed in a book entitled "Simulation for design and fabrication of VLSIs", on pp. 51-63, edited by M. Morisue and published by the CMC corporation in 1987. In this model, the following Deal-Grove equation is used. ##EQU1##
In the equation (1), t is the time, T.sub.ox is the thickness of the oxide region at the present time, T.sub.ox.sup.old is the thickness of the oxide region at the prior time, and A and B are parameters relating to the oxidation rate of a region to be oxidized.
On the other hand, individual electronic elements and/or components need to be electrically isolated in the LSI. This electrical isolation is usually realized by the selective oxidation method termed the "LOCal Oxidation of Silicon (LOCOS)" using a silicon nitride film formed on the surface of a semiconductor substrate as an oxidation mask, or the trench isolation method using trenches formed at the surface of a semiconductor substrate and filled with a dielectric.
In recent years, as the integration level of the electronic elements and components in the LSI has increased, the electronic elements and components have been miniaturized more and more. Under such the circumstances, there has been the need to simulate the isolation process for realizing the electrical isolation using the selective oxidation or trench isolation method. Also, several two-dimensional simulation methods of the isolation process have been developed.
An example of the conventional simulation methods of the isolation process using the LOCOS method is disclosed in a book entitled "Simulation Techniques of semiconductor devices and processes", on pp. 79-89, edited by K. Taniguchi and published by the Realize Incorporated in 1990. This method is explained below with reference to FIG. 1.
FIG. 1 shows the flowchart of the conventional simulation method for the LOCOS method disclosed in the Taniguchi's book.
In the step 101, desired nodes are configured onto a whole simulation zone where a SiO.sub.2 region is formed in contact with a Si region, and at the same time, the predetermined initial condition are applied to the individual nodes for setting the initial data. Also, the value of the time t is set as zero, i.e., t=0.
As seen from this description, it is assumed that the SiO.sub.2 region initially exists in contact with the Si region prior to the start of the oxidation process. In an actual oxidation process of Si, the surface of a single-crystal Si substrate is usually covered with a native SiO.sub.2 film prior to the oxidation process. Therefore, the SiO.sub.2 region is assumed to be contacted with the Si region at the start of oxidation.
In the step 102, the value of a preset time increment .DELTA.t is added to the present value (i.e., 0) of the time t. Thus, a first one of the time steps is started.
In the step 103, the following diffusion equation (2) (i.e., the Laplace's equation) of oxidant is constituted in the SiO.sub.2 region, where C.sub.ox is the concentration of the oxidant and D.sub.ox is the diffusion coefficient of the oxidant. This is because the oxidant existing in the oxidizing atmosphere is diffused through the SiO.sub.2 region to the opposing surface of the Si region. EQU D.sub.ox.cndot..gradient..sup.2 C.sub.ox =0 (2)
Then, the diffusion equation (2) is discretely solved at the individual nodes, thereby finding the surface concentration C.sub.ox.sup.surf of the oxidant at the interface between the Si and SiO.sub.2 regions (i.e., the Si/SiO.sub.2 interface) in the first time step.
Subsequently, in the step 104, using the surface concentration C.sub.ox.sup.surf of the oxidant thus found, the oxidation rate (dT.sub.ox /dt) of the Si region, which is given by the time-dependent thickness T.sub.ox of the SiO.sub.2 region, in the first time step is calculated at the individual nodes by the use of the following equation (3) ##EQU2##
where K is a coefficient of the oxidation reaction. The orientation of the oxidation rate (dT.sub.ox /dt) of the Si region is set in a direction perpendicular to the Si/SiO.sub.2 interface.
The equation (3) means that the oxidation rate of the Si region, i.e., the time-dependent thickness (dT.sub.ox /dt) of the SiO.sub.2 region, is proportional to the surface concentration C.sub.ox.sup.surf of the oxidant at the Si/SiO.sub.2 interface is assumed in this conventional simulation method.
In the step 105, a new or post-oxidation position of the Si/SiO.sub.2 interface is calculated by multiplying the value of the oxidation rate (dT.sub.ox /dt) at the Si/SiO.sub.2 interface thus found in the step 104 by the value of the time increment .DELTA.t at the individual nodes.
In the step 106, using the new or post-oxidation position of the Si/SiO.sub.2 interface thus found in the step 105, the shape or geometric deformation of the Si and SiO.sub.2 regions due to oxidation in the first time step is calculated.
In the step 107, it is judged whether the present value of the time t in the fist time step is equal to the value of the end time t.sub.END or not. If the answer is "NO", the second time step is started and the steps 102 to 106 are performed again. Further, in the same way as above, the steps 102 to 106 are repeated in the third time step and later time steps until the answer of "YES" is given. If the answer is "YES", the flow of the steps 102 to 106 is stopped.
FIGS. 2A to 2C schematically show the one-dimensional, time-dependent shape change of Si and SiO.sub.2 regions in an oxidation process, to which the above-described conventional simulation method shown in FIG. 1 is applied.
At the time t.sub.0, as shown in FIG. 2A, nodes P1, P2, P3, P4, and P5 are configured one-dimensionally along an interface L0 between Si and SiO.sub.2 regions 151 and 152. The nodes P1, P2, P3, P4, and P5 are equally spaced along the Si/SiO.sub.2 interface L0. This is performed in the step 101 in FIG. 1.
Although not shown in FIGS. 2A to 2C, it is needless to say that other nodes are configured two-dimensionally on the whole Si and SiO.sub.2 regions 151 and 152.
At the time t.sub.1 after the specific time increment .DELTA.t from the time t.sub.0 (or, in the first time step), as shown in FIG. 2B, the nodes P1, P2, P3, P4, and P5 are shifted perpendicular to the Si/SiO.sub.2 interface L0 toward the Si region 151. Thus, the nodes P1, P2, P3, P4, and P5 and the Si/SiO.sub.2 interface L0 are moved to their new positions, resulting in new nodes P1', P2', P3', P4', and P5' and a new Si/SiO.sub.2 interface L1. This movement is carried out by the use of the new or post-oxidation position of the Si/SiO.sub.2 interface L0 obtained through the steps 102 to 105 in FIG. 1.
Thus, using the result of the calculation about the oxidation rate (dT.sub.ox /dt) of the SiO.sub.2 region 152 in the step 104 and the result of calculation about the shape deformation of the Si and SiO.sub.2 regions 151 and 152 in the step 106, the thickness of the SiO.sub.2 region 152 is increased by a thickness increment .DELTA.T.sub.ox occurring in this first step.
Accordingly, the Si and SiO.sub.2 regions 151 and 152 have the shapes as shown in FIG. 2B, in which the thickness of the SiO.sub.2 region 152 is increased while the thickness of the Si region 151 is decreased due to oxidation.
At this stage, the new nodes P1', P2', P3', P4', and P5' are located on the new Si/SiO.sub.2 interface L1. The new Si/SiO.sub.2 interface L1 is apart from the old Si/SiO.sub.2 interface L0 by the thickness increment .DELTA.T.sub.ox toward the Si region 151.
This shift of the Si/SiO.sub.2 interface L0 is carried out not only when the thickness increment .DELTA.T.sub.ox is equal to or greater than a specific small value .epsilon. (i.e., .DELTA.T.sub.ox.gtoreq..epsilon.) but also when the thickness increment .DELTA.T.sub.ox is less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox &lt;.epsilon.).
Similarly, at the time t.sub.2 after the same specific time increment .DELTA.t from the time t.sub.1 (or, in the second time step), as shown in FIG. 2C, the nodes P1', P2', P3', P4', and P5' are shifted again perpendicular to the Si/SiO.sub.2 interface L1 toward the Si region 151. Thus, the nodes P1', P2', P3', P4', and P5' and the Si/SiO.sub.2 interface L1 are moved to their new positions, resulting in new nodes P1", P2", P3", P4", and P5" and a new Si/SiO.sub.2 interface L2. This movement is carried out by the use of the new or post-oxidation position of the Si/SiO.sub.2 interface L1 obtained in the steps 102 to 105.
Thus, in the same way as explained for the nodes P1, P2, P3, P4, and P5, the thickness of the SiO.sub.2 region 152 is increased by a thickness increment .DELTA.T.sub.ox ' occurring in this second time step.
Accordingly, the Si and SiO.sub.2 regions 151 and 152 have the shapes as shown in FIG. 2C, in which the thickness of the SiO.sub.2 region 152 is further increased while the thickness of the Si region 151 is further decreased due to oxidation.
At this stage, the new nodes P1", P2", P3", P4", and P5" are located on the new Si/SiO.sub.2 interface L2. The new Si/SiO.sub.2 interface L2 is apart from the old Si/SiO.sub.2 interface L1 by the thickness increment .DELTA.T.sub.ox ' toward the Si region 151.
This shift of the Si/SiO.sub.2 interface L1 is carried out not only when the thickness increment .DELTA.T.sub.ox ' is equal to or greater than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox '.gtoreq..epsilon.) but also when the thickness increment .DELTA.T.sub.ox ' is less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox '&lt;.epsilon.).
As explained above, with the above-described conventional simulation method shown in FIG. 1, to find or determine the one-dimensional, time-dependent shape change of the Si and SiO.sub.2 regions 151 and 152 shown in FIGS. 2A to 2C in each time step, the calculation steps 103 to 106 are always carried out even if the thickness increments .DELTA.T.sub.ox and .DELTA.T.sub.ox ' are less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox &lt;.epsilon. and .DELTA.T.sub.ox '&lt;.epsilon.) due to the slight thickness growth of the SiO.sub.2 region 152, as shown in FIGS. 2A to 2C.
FIG. 3 schematically shows the two-dimensional, time-dependent shape change of Si and SiO.sub.2 regions in an oxidation process, to which the above-described conventional simulation method shown in FIG. 1 is applied. FIG. 3 shows the state of the Si and SiO.sub.2 regions in the neighborhood of the edge of an oxidation mask (not shown) made of Si.sub.3 N.sub.4.
At the time t.sub.0, nodes P(n), P(n+1), P(n+2), P(n+3), and P(n+4) are configured two-dimensionally along an interface L0 between Si and SiO.sub.2 regions 141 and 143. The nodes P(n), P(n+1), P(n+2), P(n+3), and P(n+4) are spaced along the Si/SiO.sub.2 interface L0. This is performed in the step 101 in FIG. 1.
Although not shown in FIG. 3, it is needless to say that other nodes are configured two-dimensionally on the whole Si and SiO.sub.2 regions 141 and 143.
At the time t.sub.1 after a specific time increment .DELTA.t from the time t.sub.0 (or, in the first time step), the nodes P(n), P(n+1), P(n+2), P(n+3), and P(n+4) are shifted toward the Si region 141. Thus, the nodes P(n), P(n+1), P(n+2), P(n+3), and P(n+4) and the Si/SiO.sub.2 interface L0 are shifted toward the Si region 141. Thus, the nodes P(n), P(n+1), P(n+2), P(n+3), and P(n+4) and the Si/SiO.sub.2 interface L0 are moved to their new positions, resulting in new nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)' and a new Si/SiO.sub.2 interface L1. This movement is carried out by the use of the new or post-oxidation position of the Si/SiO.sub.2 interface L0 obtained through the steps 102 to 105 in FIG. 1.
Thus, using the result of the calculation about the oxidation rate (dT.sub.ox /dt) of the SiO.sub.2 region 141 in the step 104 and the result of calculation about the shape deformation of the Si and SiO.sub.2 regions 141 and 143 in the step 106, the thickness of the SiO.sub.2 region 143 is increased at the nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)' by corresponding thickness increments .DELTA.T.sub.ox occurring in this first time step.
Accordingly, the Si and SiO.sub.2 regions 141 and 143 have the shapes defined by the interface L1, as shown in FIG. 3, in which the thickness of the SiO.sub.2 region 143 is increased while the thickness of the Si region 141 is decreased due to oxidation.
At this stage, the new nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)' are located on the new Si/SiO.sub.2 interface L1. The new Si/SiO.sub.2 interface L1 is apart from the old Si/SiO.sub.2 interface L0 by corresponding thickness increments .DELTA.T.sub.ox toward the Si region 141.
This shift of the Si/SiO.sub.2 interface L0 is carried out not only when the thickness increment .DELTA.T.sub.ox is equal to or greater than a specific small value .epsilon. (i.e., .DELTA.T.sub.ox.gtoreq..epsilon.) but also when the thickness increment .DELTA.T.sub.ox is less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox &lt;.epsilon.).
Similarly, at the time t.sub.2 after the same time increment from the time t.sub.1 (or, in the second time step), the nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)' are shifted again toward the Si region 141. Thus, the nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)'and the Si/SiO.sub.2 interface L1 are moved to their new positions, resulting in new nodes P(n)", P(n+1)", P(n+2)", P(n+3)", and P(n+4)" and a new Si/SiO.sub.2 interface L2. This movement is carried out by the use of the new or post-oxidation position of the Si/SiO.sub.2 interface L1 obtained in the step 105 obtained in the steps 102 to 105.
Thus, in the same way as explained for the nodes P(n)', P(n+1)', P(n+2)', P(n+3)', and P(n+4)', the thickness of the SiO.sub.2 region 143 is increased at the individual nodes P(n)", P(n+1)", P(n+2)", P(n+3)", and P(n+4)" by corresponding increments occurring in this second time step.
Accordingly, the Si and SiO.sub.2 regions 141 and 143 have the shapes defined by the new interface L2, as shown in FIG. 3, in which the thickness of the SiO.sub.2 region 143 is further increased while the thickness of the Si region 141 is further decreased due to oxidation.
At this stage, the new nodes P(n)", P(n+1)", P(n+2)", P(n+3)", and P(n+4)" are located on the new Si/SiO.sub.2 interface L2. The new Si/SiO.sub.2 interface L2 is apart from the old Si/SiO.sub.2 interface L1 by corresponding thickness increments .DELTA.T.sub.ox toward the Si region 141.
This shift of the Si/SiO.sub.2 interface L1 is carried out not only when the thickness increment .DELTA.T.sub.ox is equal to or greater than a specific small value .epsilon. (i.e., .DELTA.T.sub.ox.gtoreq..epsilon.) but also when the thickness increment .DELTA.T.sub.ox is less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox &lt;.epsilon.).
As explained above, with the above-described conventional simulation method shown in FIG. 1, to find or determine the two-dimensional, time-dependent shape change of the Si and SiO.sub.2 regions 141 and 143 shown in FIG. 3 in each time step, the calculation steps 103 to 106 are always carried out even if the individual thickness increments .DELTA.T.sub.ox at the nodes are less than the specific small value .epsilon. (i.e., .DELTA.T.sub.ox &lt;.epsilon. and .DELTA.T.sub.ox '&lt;.epsilon.) due to the slight thickness growth of the SiO.sub.2 region 143.
With the above-described conventional simulation method as shown in FIG. 1, as explained with reference to FIGS. 2A to 2C and FIG. 3, the calculation steps 103 to 106 are always carried out independent of the magnitude of the individual thickness increments .DELTA.T.sub.ox at the individual nodes. Therefore, there is a problem that it takes very long time to simulate an oxidation process of Si.
In the case of the LOCOS method, a silicon nitride (Si.sub.3 N.sub.4) film is typically formed on a surface of a single-crystal Si substrate as an oxidation mask, and then, the Si substrate with the Si.sub.3 N.sub.4 mask is selectively oxidized in an oxidizing atmosphere. During this oxidation process, oxidant existing in the oxidizing atmosphere diffuses through a SiO.sub.2 region (i.e., native oxide of Si) existing initially on the uncovered surface of the Si substrate.
However, in the vicinity of the edge of the Si.sub.3 N.sub.4 mask, the concentration of the oxidant is very low due to existence of the Si.sub.3 N.sub.4 mask and as a result, the growth rate of SiO.sub.2 is very small.
As a result, considering this fact, it is found that unnecessary calculation is carried out in the above-described conventional simulation method shown in FIG. 1. Since the deformation calculation in the step 106 necessitates a very long time, the calculation time in the step 106 applies a large influence to the necessary simulation time.
To omit the unnecessary calculation, there is a solution that the growth rate of the SiO.sub.2 region 152 or 143 is set as zero if the thickness increment of the SiO.sub.2 region 152 or 143 is equal to or less than specific value (e.g., 1.ANG.). In this case, however, there arises the following problem.
Specifically, when the oxidation rate of the Si region 151 or 141 is low due to a low oxidizing temperature and simultaneously, the time increment .DELTA.t is set as short as possible to improve the simulation accuracy, there arises a problem that the thickness growth or increase of the SiO.sub.2 region 152 or 143 does not occur even after a long time is passed or a lot of time steps are carried out.